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ABSTRACT 



Aims. We develop a statistical analytical model that predicts the occurrence frequency distributions and parameter correlations of 
avalanches in nonlinear dissipative systems in the state of a slowly-driven self-organized criticality (SOC) system. 
Methods. This model, called the fractal-dilfusive SOC model, is based on the following four assumptions: (i) The avalanche size L 
grows as a diffusive random walk with time T, following L oc r''^; (ii) The energy dissipation rate f{t) occupies a fractal volume 
with dimension Ds, (iii) The mean fractal dimension of avalanches in Euclidean space 5 = 1,2, 3 is Z)s « (1 -I- S)/2; and (iv) The 
occurrence frequency distributions N(x) oc x^"' based on spatially uniform probabilities in a SOC system are given by N(L) oc L^^ , 
with 5 being the Eudlidean dimension. We perform cellular automaton simulations in three dimensions (S = 1,2,3) to test the 
theoretical model. 

Results. The analytical model predicts the following statistical correlations: F oc oc r^s/^ j^j. jj^g g^j,^ P oc oc r^'- for the 
peak energy dissipation rate, and E oc FT oc j^'+As/z fgj- jjjg jqj^j dissipated energy; The model predicts powerlaw distributions for all 
parameters, with the slopes Ur = (I + S)/2, af = l+(S - 1)/Ds , ffp = 2 - 1/5 , and aE = l+(S - l)KDs + 2). The cellular automaton 
simulations reproduce the predicted fractal dimensions, occurrence frequency distributions, and correlations within a satisfactory 
agreement within x 10% in all three dimensions. 

Conclusions. One profound prediction of this universal SOC model is that the energy distribution has a powerlaw slope in the range 
of Oe = 1.40 - 1.67, and the peak energy distribution has a slope of Up = 1.67 (for any fractal dimension Ds = 1, 3 in Euclidean 
space 5 = 3), and thus predicts that the bulk energy is always contained in the largest events, which rules out significant nanoflare 
heating in the case of solar flares. 

Key words. Methods: statistical - Instabilities - Sun: flares 
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1. INTRODUCTION 

The statistics of nonlinear processes in the universe often 
shows powerlaw-like distributions, most conspicously in en- 
ergetic dynamic phenomena in astrophysics (e.g., solar and 
stellar flares, pulsar glitches, auroral substorms) and in catas- 
trophic events in geophysics (e.g., earthquakes, landslides, or 
forest fires). The most widely known example is the distribu- 
tion of earthquake magnitudes, which has a powerlaw slope of 
a 2.0 for the differential frequency distribution (Turcotte 
1999), the so-called Gutenberg-Richter (1954) law. Bak, Tang, 
and Wiesenfeld (1987, 1988) introduced the theoretical con- 
cept of self-organized criticality (SOC), which has been ini- 
tially applied to sandpile avalanches at a critical angle of re- 
pose, and has been generalized to nonlinear dissipative systems 
that are driven in a critical state. Comprehensive reviews on this 
subject can be found for applications in geophysics (Turcotte 
1999), solar physics (Charbonneau et al. 2001), and astrophysics 
(Aschwanden 201 1). 

Hallmarks of SOC systems are the scale-free powerlaw dis- 
tributions of various event parameters, such as the peak energy 
dissipation rate P, the total energy E, or the time duration T 
of events. While the powerlaw shape of the distribution func- 
tion can be explained by the statistics of nonlinear processes that 
have an exponential growth phase and saturate after a random 
time interval (e.g., Willis and Yule 1922; Fermi 1949; Rosner 
and Vaiana 1978; Aschwanden et al. 1998; Aschwanden 2004, 



201 1), no general theoretical model has been developed that pre- 
dicts the numerical value of the powerlaw slope of SOC param- 
eter distributions. Simple analytical models that characterize the 
nonlinear growth phase with an exponential growth time tq and 
the random distribution of risetimes with an average value of 
ts, predict a powerlaw slope of ap - 1 + tg /tc for the energy 
dissipation rate (e.g., Rosner and Vaiana 1978; Aschwanden et 
al. 1998), but cellular automaton simulations suggest a much 
more intermittent energy release process than the idealized case 
of an avalanche with a single growth and decay phase. An alter- 
native theoretical explanation for a slope oe - 3/2 was put for- 
ward by a dimensional argument (Litvinenko 1998), which can 
be derived from the definition of the kinetic energy of convective 
flows, but this model entails a specific physical mechanism that 
has not universal validity for SOC systems. 

In this Paper we propose a more general concept where 
the powerlaw slope of the occurrence frequency distribution 
of SOC parameters depends on the fractal geometry of the 
energy dissipation domain. We aim for a universal statistical 
model of nonlinear energy dissipation processes that is indepen- 
dent of any particular physical mechanism. The fractal struc- 
ture of self-organized processes has been stressed prominently 
from beginning (Bak, Tang, & Wiesenfeld 1987, 1988; Bak 
and Chen 1989), but no quantitative theory has been put for- 
ward that links the fractal geometry to the size distribution of 
SOC events. Fractals have been studied independently (e.g., 
Mandelbrot 1977, 1983, 1985), while the fractal geometry of 
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SOC avalanches was postulated (e.g., see textbooks of Bak 
1996; Sornette 2004; Aschwanden 2011), but no general self- 
consistent model has been attempted. 

This paper presents an analytical theory that derives a the- 
oretical framework to quantitatively link the concept of fractal 
dimensions to the occurrence frequency distributions of SOC 
avalanche events (Section 2), tests of the analytical theory with 
numerical simulations of cellular automaton models in three 
Euclidean dimensions (Section 3), a comparison and application 
to solar flares (Section 4), and a summary of the model assump- 
tions and conclusions (Section 5). 

2. THEORY 

We derive in this Section a general model of the statistics of 
SOC processes, but make use of a specific example of a SOC 
avalanche that is numerically simulated with a cellular automa- 
ton code and described in more detail in Section 3, to illustrate 
and validate our theoretical derivation. 

2.1. Diffusive Random Walk in Cellular Automaton 

Avaianclies 

An avalanche in a cellular automaton model propagates via 
nearest-neighbor interactions in random directions, wherever an 
unstable node is found. The state of self-organized criticality 
ensures that the entire system is close to the instability thresh- 
old, and thus every direction of instability propagation is equally 
Ukely once a starting location is triggered (if the re-distribution 
rule is defined to be isotropic). We can therefore model the prop- 
agation of unstable nodes with a random walk in a S-dimensional 
space, which has the characteristics of a diffusion process and 
propagates in the statistical average a distance x{t) that depends 
on the time t as, 

X(0 oc _ (1) 

As a plausibility test we check this first assumption with 
an example of a numerical simulation of a cellular automaton 
model that is described in more detail in Section 3. The sim- 
ulated avalanche shown in Fig. 1 lasts for a duration of 712 
time steps and we show snapshots of the energy dissipation rate 
de/dt in Fig. 1. The complete time evolution of the energy dis- 
sipation rate de(t)/dt, the total energy e(f), the fractal dimension 
D2(t), and radius of the avalanche area r(f) as a function of time 
is shown in Fig. 3. A movie of the avalanche that shows the 
time evolution for all 712 time steps is included in the electronic 
supplementary data of this journal. The movie illustrates also 
that the instantaneous propagation direction of the avalanche is 
nearly isotropic, as we assumed here, regardless of the prior 
evolution in the neighborhood of the instantaneous energy re- 
lease. Apparendy, many of the grid points that have already been 
touched by the avalanche previously are still in a meta-stable 
state near the critical threshold that enables next-neighbor inter- 
actions. 

We calculate the total time-integrated area a{t) of the 
avalanche by summing up all unstable nodes where energy dis- 
sipation happened during the time interval [0, t] (counting each 
unstable node only once, even when the same node was unsta- 
ble more than once), and measure the mean radius r{t) of the 
avalanche area by 

r(t) ^ J—, (2) 
V n 



which closely follows the diffusive random walk distance x(t) oc 
t^^^, as it can be seen in Fig. 3 (bottom right panel), or from 
the circular area with radius r{t) drawn around the starting point 
of the avalanche in Fig. 1 (dashed circles around center marked 
with a cross). Thus, the total time-integrated area a(t) of an 
avalanche increases with a diffusive scaling. If we define T to 
be the total time duration of the avalanche, and aj = a{t = T) = 
nr^{t = T) the final area of the avalanche, then the final linear 
size L of the time-integrated avalanche is, 

L = yja^= ^r{t = T) . (3) 

According to the diffusive propagation we expect then not only 
a time evolution r{t) oc t^l^ (Eq. 1) for individual avalanches (in 
the statistical average), but also a statistical relationship between 
the size scales L and the time durations T for an ensemble of 

many avalanches, 

L oc T^l^ . (4) 

The time-integrated area ut of an avalanche, as the outUned con- 
tours in Fig. 1 show, appears to be a contiguous, space-filling 
area that is essentially non-fractal, although it has some ragged 
boundaries. For an example of a 3-dimensional avalanche see 
Mcintosh et al. (2002), which also shows essentially a space- 
filUng 3-D topology for the time-integrated avalanche volume. 
This conclusion is somewhat intuitive in cellular automaton 
models, where avalanches can propagate over nearest neighbors 
only, and thus will tend to cover a near space-filling area for 
isotropic propagation. Although the boundaries are somewhat 
ragged, the area A is equivalent to a circular area with a mean 
radius of r{t), and thus the fractal dimension would be approx- 
imately D = 2, if a box-counting method is applied within a 
quadratic area with size aj = = nr^. Thus, we define that 
the total avalanche area has a Euclidean space-filling topology 
within the diffusive boundary, and can be characterized with a 
size scale L and area aj = L^. Generalizing to 3-D space, the 
volume V of the avalanche boundary can be characterized by 
vt = L\ 

2.2. The Fractal Geometry of Instantaneous Energy 
Dissipation 

While we established the space-filling nature of the time- 
i/itegrated avalanche area with a (non-fractal) Euclidean dimen- 
sion in the foregoing section, we will now, in contrast, derive 
the theorem that the instantaneous area of energy dissipation in 
avalanches is fractal. It is actually a key concept of SOC systems 
that the spatial structure of avalanches is fractal. Bak and Chen 
(1989) express this most succintly in their abstract: "Fractals 
in nature originate from self-organized critical dynamical pro- 
cesses". 

As it can be seen from the snapshots of an evolving 
avalanche shown in Fig. 1, the instantaneous areas of energy 
dissipation cover a fraction of the solid area a{t) that is encom- 
passed by the diffusive boundary. A detailed inspection of the 
shapshots shown in Fig. 1 even reveals a checkerboard pattem of 
instantaneous avalanche maps that emphasizes the fractal topol- 
ogy of cellular automaton avalanches. We make now the sec- 
ond major assumption that the area A(t) of instantaneous energy 
dissipation de{t)ldt is fractal, or that the volume V{t) for 3-D 
avalanches is fractal, respectively. To simplify the nomenclature, 
we will generally refer to the fractal volume Vs in S-dimensional 
Euclidean space, which corresponds to V'3 = V for fractal vol- 
umes in 3-D space, to Vi = A for fractal areas in 2-D space. 
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Fig. 1. Time evolution of the largest avalanche event #1628 in the 2-D cellular automaton simulation with grid size - 64^. The 
12 panels show snapshots at particular burst times from f = 26 to f = 690 when the energy dissipation rate peaked. Active nodes 
where energy dissipation occurs at time t are visualized with black and grey points, depending on the energy dissipation level. The 
starting point of the avalanche occurred at pixel ix,y) = (41,4), which is marked with a cross. The time-integrated envelop of the 
avalanche is indicated with a solid contour, and the diffusive avalanche radius r{t) = t^^^ is indicated with a dashed circle. The 
temporal evolution is shown in a movie available in the on-line version. 
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Fig. 2. Determination of the fractal dimension D2 - 
log A, / log X, for the instantaneous avalanche sizes of the 12 time 
steps of the avalanche event shown in Fig. 1. Each row is a dif- 
ferent time step and each column represents a different binning 
of macropixels (Ax, = 1,2,4,8). The fractal dimension is de- 
termined by a linear regression fit shown on the right-hand side. 
The mean fractal dimension of the 12 avalanche snapshots is 
D2 = 1.43 ±0.17. 



and Vi = X for fractal lengths in 1-D space. (Note that we use 
uppercase symbols V,A,X for fractal parameters, while we use 
lowercase symbols v, a, x for non-fractal Euchdean parameters). 



A fractal volume Vs can be defined by the HausdorfF dimension 
Ds in S-dimensional Euclidean space, 



Ds - lim 



logVj 



^0 logx 



(5) 



where Vs is the fractal volume with Euclidean scale x, or by the 
scaling law. 



(6) 



In Fig. 2 we demonstrate the fractal nature of the 12 instan- 
taneous avalanche snapshots shown in Fig. 1. We rebin the 
avalanche area into macropixels with sizes of Ax, = 2',i = 
0, ...,3 (or X = 1,2,4,8), measure the number of macropixels 
Ai that cover the instantaneous avalanche area, use the diffu- 
sive scaling x(0 oc (Eq. 1) for the rebinned length units 
X,- = x{t)/Axi, and determine the Hausdorff dimemsion D2 from 
the linear regression fit log(A,) - D2 log(x,). We find that that the 
4 datapoints for each of the 12 cases exhibit a linear relationship, 
which proves the fractaUty of the avalanche areas. The average 
fractal dimension of the 12 timesteps shown in Fig. 1 and 2 is 
D2 = 1.43 + 0.17. 

We measure now the fractal dimension D2 of the instanta- 
neous energy dissipation volume in the 2-D avalanche for all 700 
time steps of its duration, shown for 12 time instants in Fig. 1. 
The unstable nodes (signifying instantaneous energy dissipation) 
are counted in each time step, which yield a number for the frac- 
tal volume or area A(t) = ¥2(1), while the size x of the encom- 
passing box is determined from the area of the time-integrated 
avalanche, i.e., x{t) = ^/a(fj, which yields the time evolution of 
the fractal dimension D2{t) = log[V2(f)]/log[x(f)] (Eq. 5) as a 
function of time, shown in the top right panel in Fig. 3. The frac- 
tal dimension £>2(0 fluctuates around a constant mean value of 
D2 - 1.45 + 0.13, which is close to the arithmetic mean of the 
minimum dimension D2.,nin ~ 1 and maximum Euclidean limit 
D2,max = 2, i.e. {D2} ~ (£>2,mH7 + £>2,mfl.t)/2 = 3/2. TMs cor- 
roborates our second major assumption that the instantaneous 
volume of energy dissipation is fractal, and that the fractal di- 
mension can be approximated by a mean (time-independent and 
size-independent) constant during the evolution of avalanches, 
in the statistical average. 

If we moreover define a mean energy dissipation rate quan- 
tum (AE) per unstable node, which is indeed almost a constant 
for a cellular automaton model near the critical state, we expect 
a scaling of the instantaneous dissipation rate (or flux) f(t) that 
is proportional to the instantaneous dissipation volume Vs (with 
Eq. 6), 



f(t) = ^ cx (AE) Vs(t) = (AE) x(tf^ 
at 



(7) 



Combining this with the diffusive expansion of the boundary 
x{t) oc t^/"^ (Eq. 1), we can then predict the average time evo- 
lution of the energy dissipation rate fit) = de{t)ldt. 



/(0 = ^cx(A£)/^^/^) 
at 



(8) 



Integrating Eq. (8) in time, we obtain the time evolution of 
the total dissipated energy, e{t). 



Jo dT Jo 



>(l+£>s/2) 



(9) 



Hence, for our 2-D avalanche (with D2 - 3/2) we expect an evo- 
lution of e{t) oc which indeed closely matches the actually 
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Fig. 3. Time evolution of the largest avalanche event #1628 in the 2-D cellular automaton simulation with grid size - 64^. The 
time profiles include the instantaneous energy dissipation rate f(t) = de/dt (top left), the time-integrated total energy e(t) (bottom 
left), the instantaneous fractal dimension D2(t) (top right), and the radius of the avalanche area r{t) (bottom right). The observed 
time profiles from the simulations are outhned in solid hnestyle and the theoretically predicted average evolution in dashed hnestyle. 
The statistically predicted values of the instantaneous energy dissipation rate fit) cc t^l^ (dotted curve) and peak energy dissipation 
rate p{t) oc (dashed curve) after a time interval t are also shown (top left panel). The 12 time labels from 26 to 690 (top left frame) 
correspond to the snapshot times shown in Fig. 1. 



simulated cellular automaton case, as we see in Fig. 3 (bottom 
left panel). The time evolution of the energy dissipation rate is 
shown in Fig. 3 (top left panel), which fluctuates strongly dur- 
ing the entire avalanche, but follows in the statistical average 
the predicted evolution de(t)/dt oc f^^/^ ^ p/4) f^j- ^ 3/2. 
Note that our analytical expressions of the time evolution of 
avalanches, such as the hnear size x(t) (Eq. 1), the instantaneous 
energy dissipation rate f(t) (Eq. 7, 8), or the total dissipated en- 
ergy e(t) (Eq. 9), do not predict the specific evolution of a single 
avalanche event, but rather the statistical expectation value of a 
large ensemble of avalanches, similar to the statistical nature of 
the difiusive random walk model (Eq. 1). 

The time evolution of the instantaneous energy dissipation 
rate fit) fluctuates strongly, as it can be seen in for the largest 



avalanche simulated in a cellular automaton model (Fig. 3, top 
left panel). We might estimate the peak values that can be ob- 
tained statistically (after a time duration t) from the optimum 
conditions when the fractal filling factor of the avalanche reaches 
a near-Euclidean filling, i.e., in the limit of Ds i-> 5. Replacing 
the fractal dimension Ds by the Euclidean limit S in Eqs. 7 and 8 
yields then (as an upper hmit) an expectation value for the peak 
Pit), 

pit) oc <A£) V^''\t) oc <A£) xit)^ oc <A£) t*-^'^^ . (10) 

Denoting the energy dissipation rate in the statistical aver- 
age after time t = T with with F = fit = T), the peak energy 
dissipation rate with P = pit - T), and the total energy of the 
avalanche with E = eit = T), it follows from Eqs. (8-10) that 
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E oc FT - P^' 1^ T, and we expect then the following correla- 
tions between the three parameters E, F, P and T for an ensem- 
ble of avalanches, 



F oc 
p oc 



(11) 



For instance, for a 2-D avalanche with an average fractal dimen- 
sion of D2 = 3/2 we expect the following two correlations, 
E oc T''^*, F oc T^^'^, and P oc r' (see Fig. 3). The power- 
law indices for the correlated parameters are listed for the three 
Euclidean dimensions S = 1, 2, 3 separately in Table 1. 

2.3. Occurrence Frequency Distributions 

Considering the probability of an avalanche with volume V, the 
statistical likelihood simply scales reciprocally to the volume 
size V, if avalanches are equally likely in every space location 
of a uniform volume Vq of a system in a (self-organized) critical 
state. This is illustrated in Fig. (4). For the 1-D Euclidean space, 
n = I avalanche can happen with the maximum size L = Lq 
of the system (top left), n = 2 avalanches with the half size 
L - Lo/2, or n = 4 avalanches for a quarter size L = Lq/4. For 
the 2-D Euclidean space, the number of possible avalanches that 
can be fit into the total area Aq of the system is n = 1 for A = Aq, 
„ = 2^ = 4 for A = Ao/2, or n = 2"* = 16 for A = Ao/4 (second 
row). For the 3-D Euclidean space we have, correspondingly, 
n = I for y = Vo, n = 2^ = 8 for cubes of half size L = Lo/2, 
and n = 4^ = 64 for quarter-size cubes with L = Lo/4 (bottom 
row). So, generalizing toS = 1, 2, 3 dimensions, we can express 
the probability for an avalanche of size L and volume Vs = 
as. 



NiL) oc ' oc l '^ 



(12) 



This simple probability argument is based on the assumption 
that the number or occurrence frequency of avalanches is equally 
likely throughout the system, so it assumes a homogeneous dis- 
tribution of critical states across the entire system. 

The occurrence frequency distribution of length scales, 
N{L) oc (Eq. 12), serves as a primary distribution function 
from which all other occurrence frequency distribution functions 
A^(x) can the derived that have a functional relationship to the 
primary parameter L. First we can calculate the occurrence fre- 
quency distribution of avalanche time scales, by using the diSu- 
sive boundary propagation relationship L{T) oc r'^^ (Eq. 4), by 
substituting the variable T for L in the distribution N(L) (Eq. 1 2), 



NiT)dT = NiL[T]) 



dL 
df 



dT oc r-«i^^)/2] ^ 



(13) 



Subsequently we can derive the occurrence frequency distribu- 
tion function N{F) for the statistically average energy dissipa- 
tion rate F = f{t = T) using the relationship F{T) oc T^s/^ 
(Eq. 11), 



N{F)dF = N{T[F]) 



dT 



dF 



(14) 



the occurrence frequency distribution function of the peak en- 
ergy dissipation rate P using the relationship relationship P{T) oc 
T^l^ (Eq. 11), 



N{P)dP = N{T[P]) 



dT 



dP 



dP oc p-P-i/s] ^ 



and the occurrence frequency distribution function N{E) for the 
total energy E using the relationship E{T) oc (Eq. 1 1), 



N{E)dE = NiT[E]) 



dT 



dE 



dE oc £-[i+(5-i)/(i5s+2)] _ 



(16) 



Interestingly, this derivation yields naturally powerlaw functions 
for all parameters L, T, F, P, and E, which are the hallmarks 
of SOC systems. In summary, if we denote the occurrence fre- 
quency distributions N(x) of a parameter x with a powerlaw dis- 
tribution with power index Ox, 



N{x)dx oc X dx , 



(17) 



we have the following powerlaw coefiicients Ux for the parame- 
ters X = T,F, P, and E, 



ax = (l+S)/2 

ap - 1 + (5 - l)/Ds 

ap = 2- l/S 

ap =1+{S- l)KDs + 2) 



(18) 



For instance, for our 2-D cellular automaton model with 5=2 

and Ds -3/2 we predict powerlaw slopes of aj- = 3/2 = 1.5, 
ap = 5/3 ~ 1.67, ap = 3/2 = 1.5, and ap = 9/7 ^ 1.28. The 
powerlaw coefficients are summarized in Table I separately 
for each Euclidean dimension 5 = 1,2,3. 

2.4. Estimating the Fractal Dimension of Cellular Automatons 

The dynamics of SOC systems is often simulated on computers 
with cellular automaton codes, which have a S-dimensional lat- 
tice of nodes, where a discretized mathematical redistribution 
rule is applied once a local instability threshold is surpassed 
(e.g., Bak et al. 1987; Lu and Hamilton 199 1; Charboimeau et 
al. 2001). Avalanches in such cellular automaton models propa- 
gate via nearest-neighbor interactions, which includes (25 + 1) 
nodes in a S-dimensional lattice grid: one element is the unsta- 
ble node, and there are 2S next neighbors. For instance, a 1- 
dimensional grid in Euclidian space with dimension 5 = 1 has 
(25 -H 1) = 3 nodes involved in a single nearest-neighbor relax- 
ation step, a 2-dimensional grid (5 = 2) has (25 -H 1) = 5 nodes, 
and a 3-dimensional grid (5 = 3) has (25 -n 1) = 7 nodes (Fig. 5). 

How can we estimate the fractal dimension for avalanches 
in such a lattice-based cellular automaton model? The minimum 
fractal dimension of a growing 3-D avalanche corresponds to a 
1-D Unear structure, D3 ,,,,,, = 1, because avalanches evolve by 
iterative propagation from one node to the next-neighbor node, 
so a contiguous linear path from one node to the next neighbor 
is about the sparsest spatial structure that still enables avalanche 
growth in a SOC model. If the linear avalanche path would be 
discontinuous (Ds < 1), the avalanche is likely to die out due to 
a lack of unstable next neighbors, so Ds^i„ ~ 1 represents prac- 
tically a lower cutoff", which does not exclude occasional values 
of Ds < 1 (fractal dust). At the other extreme, when all nodes 
are close to the instability threshold, an avalanche can aff"ect all 
next neighbors and grow as a nearly space-filling structure with 
a maximum fractal dimension of Ds,max = S that equals the 
Euclidean space dimension 5 . The mean average fractal dimen- 
sion {Ds ) can then be estimated from the geometric mean of the 
minimum Vs^min and maximum fractal volume Vs^max, which is 
equivalent to the arithmetic mean of the minimum Ds^min and 
maximum fractal dimension Ds^x (Eq. 5), 



(15) (Ds) 



l0g{V) log V^^^i^Jin^^^^ Ds,mm + D, 



S,max 



logx 



logx 



(19) 
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Table 1. Theoretically predicted occurrence frequency distribution powerlaw slopes a and power indices jS of parameter correlations 
predicted for SOC cellular automatons with EucUdean space dimensions S = 1, 2, 3. 



Parameter 


Theory 


5 = 


1 


S = 2 


5=3 




Fractal Dimension: 


Ds =(l+5)/2 






3/2 


2 


Length scale powerlaw slope: 


tti = S 






2 


3 


Duration powerlaw slope: 


ar = (1 +5)/2 






3/2 


2 


Instantaneous energy dissipation rate slope: 


aF = l+(S- l)/Ds 






5/3 


2 


Peak energy dissipation rate slope: 


ap = 2- IIS 






3/2 


5/3 


Energy powerlaw slope: 


aE = \+{S- l)l{Ds + 2) 






9/7 


3/2 


Diffusive scaling of length L with duration T, 


L oc 


L oc 


jl/2 


Locri/2 


L oc 


Correlation of peak rate F with duration T , 


F oc 


Foe 


jl/2 


F oc T^!* 


Foe r' 


Correlation of peak rate P with duration T , 


p oc 


P oc 


7^1/2 


Poc r' 


p oc r3/2 


Correlation of energy E with duration T , 


E oc ji+Bs/z 


E oc 


J'3/2 


E oc r7/4 


EocT^ 



From this we predict the following mean fractal dimensions 
{Ds) in different Euclidean spaces with dimensions S = 1, 2, 3, 



<Di> = (l + l)/2= 1 
(D2) - (1 + 2)/2 = 3/2 
(£)3> = (l+3)/2 = 2. 



(20) 



or more generally as a function of the Euclidean space dimension 
S, 



(Ds) 



forS = 1,2,3 , 



(21) 



Using these estimates of the mean fractal dimension {Ds ) into 
Eq. (18) we obtain the numerical values given in Table 1. Note, 
that this estimate of the mean fractal dimension is only an ap- 
proximation, because the lower limit „,,„ « 1 is not rigor- 
ously derived from probability theory. However, as the example 
in Fig. 3 shows, the estimated fractal dimension of D2 = 1.5 
comes close to the observed mean value of {D2} = 1.45 ± 0.13 
in this avalanche. 

SOC systems that are different from the isotropic cellular 
automaton model we are using here may have different fractal 
dimensions. Thus our theoretical prediction of the mean fractal 
dimension applies only to SOC systems with similar isotropic 
next-neighbor redistribution rules. For any observed SOC sys- 
tem, the fractal dimension Ds can be empirically determined by 
measuring the powerlaw slopes af or Oe, which are a function 
of the fractal dimension Ds (Eq. 18). 

3. CELLULAR AUTOMATON SIMULATIONS 

3.1. Numerical Cellular Automaton Code 

An isotropic cellular automaton model that miinics a SOC sys- 
tem was originally conceived by Bak et al. (1987, 1988) and first 
applied to solar flares by Lu and Hamilton (1991). A version 
generalized to 5 = 1,2,3 dimensions is given in Charbonneau 
et al. (2001) and is also summarized in Aschwanden (201 1). 

The numerical simulations of cellular automaton models 
take place in a S-dimensional cartesian grid, where nodes are 
localized by discretized coordinates, say with Xijk in S =3, and 
a physical scalar quantity Bij^ = B(x = Xijk) is assigned to each 
node. The dynamics of the system is initiated by quantized en- 
ergy inputs 6B{t) at random locations Xijk with a constant rate as 
a function of time. At each time step t and spatial node xijk, the 
local stabiUty is checked by measuring the local S-dimensional 



"curvature" (Charbonneau et al. 2001) with respect to the next 
neighbor cells (nodes) x„„ = Xi±ij±i^k±u 



AS, 



ijk 



B, 



ijk ' 



2S Tj^""- 



(22) 



Most nodes are stable at a particular time if the system is driven 
slowly, say with an input rate 6BI{Bijk) 1 (Charbonneau et 
al. 2001). The system is defined to be stable, as long as the lo- 
cal gradients are smaller than some critical threshold value Be. 
However, once a local gradient exceeds the critical value, i.e., 
^ijk ^ Be, a mathematical redistribution rule is applied that 
smoothes out the local gradient and makes it stable again. The 
redistribution rule simply spreads the difference ABij^ (or the 
threshold B^) equally to the next neighbors (in an isotropic cel- 
lular automaton model). 



2iS 1-» T» 

Note that the amount of the redistributed quantity is the thresh- 
old energy Be in the models of Lu et al. (1993) and Charbonneau 
et al. (2001), while it is the (larger) amount of the actual gradi- 
ent ABijk in the original model of Lu & Hamilton (1991). If the 
node is unstable, then the actual gradient is larger than the crit- 
ical gradient, rather than smaller. This was modified in a later 
paper (Lu et al. 1993) by redistributing the threshold gradient 
rather than the full gradient, presumably due to numerical insta- 
bilities (Liu et al. 2002). This redistribution rule is conservative, 
in the sense that the quantity B is conserved after every redistri- 
bution step, because the same amount is transferred to the next 
neighbors that is taken away from the central cell. However, al- 
though the scalar field quantity B is conserved, the energy B^ is 
not conserved after a redistribution step, because of the nonlin- 
ear (quadratic) dependence assumed, which was introduced in 
analogy to the magnetic field energy density Emag = B^ISn. In 
fact, every redisttibution of \AB\ > Be dissipates energy from the 
system, by an amount of (Charboimeau et al. 2001), 



-^Be. 

25-1-1 



(23) 



(24) 



Thus, the minimum amount of dissipated energy is for |AS| > Be, 
when the threshold Be is infinitesimally exceeded by \AB\, 

Emin = ^Hy^' • (2^) 

Once a node Xjjk is found to be unstable, the check of unsta- 
ble cells propagates to the next neighbors x^n — Xi±ij±i,k±i 
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1-D Avalanches 



I r 



x=1 
n=1 
L=1 



x=1/2 
n=2 

L=x^=(1/2)^ 



2-D Avalanches 




x=1/2 

n=2^=4 

A=x'=(1/2)' 



x=1/4 
n=4 

L=x^=(1/4)^ 



x=1/4 

n=4^=16 

A=x2=(1/4)' 



3-D Avalanches 




x=1 
n=1 
V=1 



x=1/2 

n=2^=8 

V=x^=(1/2)^ 



x=1/4 

n=4^=64 

V=x^=(1/4)^ 



Fig. 4. Schematic diagram of the Euclidean volume scaling of the diffusive avalanche boundaries, visualized as circles or spheres 
in the three Euclidean space dimensions S - 1, 2, 3. The Euclidean length scale x of subcubes decreases by a factor 2 in each step 
(x, = 2 ' , / = 0, 1 , 2), while the number of subcubes increases by n, = (2')'^ , defining a probability of A^(x,) oc for each avalanche 
size with size x, . 



the next time step and all unstable neighbor cells are subject to 
the redistribution rule, and progressively continues to the next 
neighbors each time step until all cells are stable again. Such 
a chain reaction of next-neighbor redistributions is called an 
avalanche event (see examples in Fig. 4). Note that a minimum 
avalanche has to include at least one redistribution step. 

In this study the author coded independently such a cellu- 
lar automaton algorithm according to the specifications given in 
Section 2.5 of Charbonneau et al. (2001) and it was run with 
exactly the same system parameters, such as the input quan- 
tity cr\ < 6B < 0-2 homogeneously distributed in the range of 
(Ti = -0.2 to 0-2 - 0.8, the threshold quantity Be - 5, with grid 
sizes of N = 128 and 256 in one dimension (S = 1), = 32 
and 64 in two dimensions (S = 2), and A' = 16 and 24 in three 



dimensions (5 = 3). We sampled the total volumes V, energies 
E, peak energy dissipation rate P, and durations T of avalanches 
and were able to reproduce the results given in Charbonneau et 
al. (2001) consistently, although we used different random gen- 
erators for the input, different time intervals for the initiation 
phase and onset of SOC (tsoc = 3 x 10^ for = 128 and 5 = 1; 
tsoc = 15 X 10^ for A? = 256 and 5 = 1; tsoc = 5 x 10*^ for 
A^ = 32 and 5 = 2; tsoc = 41 x lO*" for A^ = 64 and 5 = 2; 
tsoc =4x10*" for A' = 16 and 5 = 3; tsoc = 15 x 10^ for 
N = 24 and S =3), and slightly different powerlaw fitting pro- 
cedures. In addition, we calculated also the fractal dimensions of 
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S=1 
nn=3 





Fig. 5. The fractal geometry of next-neighbor interactions are shown for a cellular automaton lattice model for Euclidian dimensions 

S - I (left), 5=2 (middle), and 5=3 (right). The unstable node is shaded with grey, and the next neighbor nodes in white. The 
total number nn of nodes involved in a local redistribution rule scales as nn = 2S + I with the EucUdean dimension S . 



the SOC avalanches from the avalanche volumes V and the size 
X of the enveloping EucUdean cube , 



log(V) 
log(x) ' 



(26) 



where x is the largest spatial scale that brackets the fractal 
avalanche volume in each spatial direction (x,y,z) of the S- 
dimensional Euclidean space. This definition is slightly different 
from the "radius of gyration" method employed in Charbonneau 
et al. (2001), but follows more the standard convention of frac- 
tal dimensions measured with box-counting methods. In the fol- 
lowing we show the results from two runs in each dimension 
S = 1,2,3 and compare them with the theoretical predictions 
made in Section 2. 



3.2. Occurrence Frequency Distributions 

The results of occurrence frequency distributions and correla- 
tions are shown for a 1-D (Fig. 6), a 2-D (Fig. 7), and a 3-D 
cellular automaton code (Fig. 8), and listed in Table 2. We eval- 
uated the powerlaw slopes by a weighted linear regression fit 
(weighted by the number of avalanche events per bin) by ex- 
cluding undersampled bins (less than 20 events). Since the exact 
values depend sometimes on the fitted range, we vary the range 
of linear regression fits from the full range of the powerlaw part 
(indicated with black line in Figs. 6-8) to the upper half (indi- 
cated with grey line in Figs. 6-8) and calculate the average and 
standard deviations of the powerlaw slope values. 

The values of the predicted and simulated powerlaw slopes 
are compiled in Table 2. First of all we notice an extremely good 
agreement (within a few percents) of the powerlaw slopes ue, 
ap, and aj- between our simulations and those of Charbonneau 
et al. (2001) and Mcintosh et al. (2002), although we used dif- 
ferent codes, initiation times tsoc, and powerlaw fitting pro- 
cedures. We performed the 1-D to 3-D runs with small cubes 
(A^ = 128, 32, 16) as well as with larger cubes (A^ = 256, 64, 24), 
which both give very consistent values (Table 2). In the follow- 
ing we quote only the values for the larger cubes, which are also 
shown in Figs. 6-8. 

Considering the agreement between theory and simulations, 
there is a reasonable agreement for all three space dimensions 
S = 1,2,3. For the fractal dimension of avalanches we find 



D2 = 1.60+0. 17 (predicted D2 = 1.5)and/)3 = 1.94 ±0.27 (pre- 
dicted D3 = 2.0), which are fully consistent with our theoretical 
estimate of Ds = (1 +S)/2 (Eq. 21). 

The relationship of the avalanche probability being recipro- 
cal to the volume, N(L) oc (Eq. 12) is also approximately 
confirmed by the simulations, for which we find ai = 0.88±0.09 
(predicted = 1 for 5 = 1), ai = 2.14 ± 0.18 (predicted 
az, = 2 for 5 = 2), and - 2.55 ± 0.10 (predicted = 3 
for S = 3), where the latter value has the largest error due to the 
smallest sizes of 3-D cubes with a dynamic range of only about 
1 dex. Clearly this relationship agrees better for larger cubes (up 
to L = 256 in 1-D simulations). We have to add a caveat that 
finite-size effects are likely to restrict the size of avalanches at 
the system boundaries, especially for the 3D cellular automata 
runs which we run here with a size of 16^ and 24^ only. Also un- 
certainty estimates of the powerlaw slopes are less accurate for 
small dynamic ranges (i.e., small system size L here) and due to 
histogram binning. In the real world, however, finite-size effects 
may also cause modifications of powerlaw distributions, such as 
the maximum size of active regions or the vertical density scale 
height of coronal loops. 

The total energy powerlaw slope qe (Eq. 18) fitted from the 
powerlaw distributions A^(£') yields values of ue = 1.06 ± 0.04 
in 1-D (predicted oe = 1.0), oe = 1.48 ± 0.03 in 2-D (predicted 

= 1.28), and = 1.50 ± 0.06 in 3-D (predicted Qf£ = 1.5), 
which agree with the predictions within 6%, 16%, and 0%. 

The peak energy dissipation rate powerlaw slope ap (Eq. 18) 
fitted from the powerlaw distributions N(P) yields values of 
ap = 0.94 ± 0.18 in 1-D (predicted ap = 1.0), ap = 1.85 ± 0.06 
in 2-D (predicted ap = 1.5), and ap = 1.96 + 0.14 in 3-D (pre- 
dicted ap = 1.67), which agree with the predictions within 6%, 
23%, and 17%. 

The time duration powerlaw slope ar (Eq. 18) fitted from the 
powerlaw distributions Nil) yields values of = 1.17 ± 0.02 
in 1-D (predicted Ot = 1.0), ar = 1.77 ± 0.18 in 2-D (predicted 
ap = 1.5), and ap = 1.76 + 0.19 in 3-D (predicted ap = 2.0), 
which agree with the predictions within 17%, 18%, and 12%. 

Most occurrence frequency distributions are subject to a 
drop-off from an ideal powerlaw distribution due to finite-size 
effects. It is quite satisfactory that our simple first-order theory 
predicts most powerlaw slopes within an accuracy of order 10%. 
We have also to be aware that our first-order theory assumes that 
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S0C1 _Q0_256_01 5jDar.dat 
ttp = 0.94+0.18 



N= 62968 
aL= 0.88+0.09 



10 
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10 

Length scale L 




(iflfSXT.*-. :.• 




100 



|l!'lJif^70:67+0.20 
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Total duration T(s) 



P*T= E 



0.99+0.02 



10 100 
Total duration T(s) 
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10^ 
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1000 
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1000 



10 

Length scale L 



100 



Fig. 6. Cellular automaton simulations with a N = 256 1-D lattice produced by a numerical code according to Charbonneau et 
al. (2001). The frequency distributions of peak energy dissipation rate P, total energies E, time durations T, and fractal avalanche 
volumes V are shown along with the fitted powerlaw slopes. Correlations between the fractal dimension D\ and parameters E, P, 
and T are also shown, fitted in the ranges of P > 5Q, E > 50, T > 5, and V >2. Only a representative subset of 2000 events are 
plotted in the scatterplots. 



the fractal dimension of the energy release rate is a constant, 
while the simulated avalanches show a slight trend of increasing 
fractal dimensions D with size L, i.e., D2 ~ LP'^ for 5=2, 
and D3 X L^-^^ for 5=3 (Figs. 7 and 8, bottom right). This 
slight trend of Dg (L) affects the powerlaw slopes ap and to 
be somewhat flatter for small avalanches than for larger ones, 
which represents a second-order effect and could be considered 
in a more refined theory. 



3.3. SOC Parameter Correlations 

An alternative method of testing our theory is (i) a determina- 
tion of the power index yS of correlated parameters x and y, i.e., 
>' oc ^c^, by linear regression fits log(j) oc log(>'o) -1- /S log(x), or (ii) 
by inferring them from the powerlaw slopes and ay of their 
occurrence frequency distributions, i.e., /3 - (a^ - l)/(ay - 1) 
(see derivation in Section 7.1.6 of Aschwanden 2011). The re- 
sulting values are listed for both methods in Table 3 (labeled 
with "Regression" and "Slopes"), for each 1-D, 2-D, and 3-D 
simulation run of our cellular automaton code. The linear re- 
gression fits are shown in the lower halves of the Figs. 6, 7, and 
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Fig. 7. Cellular automaton simulations with a N - 64^ 2-D lattice produced by a numerical code according to Charbonneau et 
al. (2001). Representation otherwise similar to Fig. 6. 



8. Note that truncations occurs for each parameter due to the ef- 
fect of finite-size systems that have been used in the numerical 
simulations (L =128 and 256 for 1-D; L =32 and 64 for 2-D; and 
L -16 and 24 for 3-D lattices). The minimum amount of dissi- 
pated energy for a threshold of fi^ = 5is£TO„ = 16.7 (for 5 = 1), 
E^in = 20 (for S = 2), and E„i„ = 21 A (for S = 3), according 
to Eq. (26). Thus we perform the linear regression fits only in 
parameter ranges that are not too strongly aff'ected by truncation 
effects, which is at durations T > 5, peak energy dissipation rate 
P > 50, and energies E > 50. To quantify an error of the Unear 
regression fits, we perform fits of y{x), and with exchanged axes, 
x(y), and quote the mean and half difl'erence for the two fits. 

The correlation L oc r'^^ (Eq. 4) tests our assumption of 
a diffusive random walk for the propagating avalanche bound- 



aries. For the theoretically expected value of the power index 
Ptl = 0.5 we find ^tl = 0.53 - 0.65 for 1-D avalanches, 
Ptl = 0.61 - 0.76 for 2-D avalanches, and fin = 0.48 - 0.53 
for 3-D avalanches, which corroborates our assumption of a dif- 
fusive avalanche expansion. 

For the correlation of the peak energy dissipation rate with 
the duration of an avalanche P oc T^^p fjj^^j - q 57 _ 0.73 
for 1-D avalanches (predicted ySr^ = 0.5), /Stp = 0.91 - 1.04 for 
2-D avalanches (predicted /3tp = 1-0), and /3tp = 0.79 - 1.02 
for 3-D avalanches (predicted firp = 1-5), which amounts to an 
agreement of « 30% for the 1-D case, « 3% for the 2-D case, and 
a: 33% for the 3-D case. We suspect that a lot of avalanches are 
stopped at the boundary, which underestimates the peak energy 
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Peak counts P{cts/s) Length scale L Total duration T(s) 




Fig. 8. Cellular automaton simulations with a N - 24^ 3-D lattice produced by a numerical code according to Charbonneau et 
al. (2001). Representation otherwise similar to Figs. 6 and 7. 



dissipation rate and thus yields a systematically too high power 
index fijp for the 3-D case. 

For the correlation of the total energy with the duration of 

an avalanche E oc T^'te we find /Ste - 1.56 - 1.65 for 1-D 
avalanches (predicted /3rp = 1.5), fijE - 1-61 - 1.82 for 2- 
D avalanches (predicted Ptp - 1-75), and Pte - 1.51-1.75 
for 3-D avalanches (predicted Ptp = 2.0), which amounts to an 
agreement of « 1% for the 1-D case, « 2% for the 2-D case, and 
« 18% for the 3-D case. 

In summary, we find an overall agreement of order 10% be- 
tween theory and numerical simulations for the power indexes 
of correlated parameters. 



4. OBSERVATIONS OF SOLAR FLARES 

Let us turn now to some astrophysical observations to compare 
our theory and the results of cellular automaton simulations. Lu 
and Hamilton (1991) applied SOC theory and cellular automa- 
ton models for the first time to solar flares. They modeled the 
flare statistics from hard X-ray counts observed with the HXRBS 
detectors onboard SMM, which shows a powerlaw distribution 
extending over 4 orders of magnitude, with a powerlaw slope 
of ap =s 1.8 for the peak count rate P (Dennis 1985). From 
their cellular automaton code they obtained powerlaw slopes of 
ap » 1.8 for the peak energy dissipation rate P, aj ~ 1.9 for 
the flare durations T, and oe ~ lA for the total energies E. We 
re-analyzed the HXRBS/SMM dataset and obtain the values of 
ap * 1.73 ± 0.01 for the power P, aj ~ 2.29 ± 0.25 for the 
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Table 2. Theoretically predicted and numerically simulated powerlaw slopes of occurrence frequency distributions of cellular au- 
tomaton models with Euclidean dimension S = 1, 2, 3 in the state of self-organized criticality. 



Reference 


Theory 


S = 1 


5=2 


5=3 





Fractal Dimension: Ds 



Theory: Z)s=(l+S)/2 


1.00 


1.50 


2.00 


Simulation (N=256,64,24) 


1.00 ±0.00 


1.60 ±0.17 


1.94 ±0.27 


Simulation (N=128, 32,16) 


1.00 ±0.00 


1.62 ± 0.18 


1.97 ± 0.29 


Charbonneau et al. (2001) 




1.58 ± 0.02 


1.78 ± 0.01 


Mcintosh et al. (2002) 




1.58 ± 0.03 


1.78 ± 0.01 


Length scale powerlaw slope: ol 








Theory: ai=S 


1.00 


2.00 


3.00 


Simulation (N=256,64,24) 


0.88 ± 0.09 


2.14 ±0.18 


2.55 ± 0.10 


Simulation (N=128, 32,16) 


0.89 ± 0.10 


1.96 ± 0.09 


2.58 ± 0.11 


Energy powerlaw slope: Ue 








Theory: q-£ 


1.00 


1.28 


1.50 


Simulation (N=256,64,24) 


1.06/ ±0.04 


1.48 ± 0.03 


1.50 ± 0.06 


Simulation (N=126,32,16) 


1.09 ± 0.05 


1.4/ ± U.U3 


1.51 ± 0.06 


Charbonneau et al. (2001) 




1.42 ± 0.01 


1.47 ± 0.02 


Mcintosh ct al. (2002) 




1.41 ± 0.01 


1.46 ± 0.01 


Peak energy rale powerlaw slope: ap 








Theory: ap 


1.00 


1.50 


1.67 


Simulation (N=256,64,24) 


0.94 ±0.18 


1.85 ±0.06 


1.96 ± 0.14 


Simulation (N=128,32,16) 


1.05 ±0.10 


1.73 ± 0.08 


1.95 ±0.13 


Charbonneau et al. (2001) 




1.72 ±0.02 


1.90 ±0.03 


Duration powerlaw slope: ot 








Theory: ar 


1.00 


1.5 


2.00 


Simulation (N=256,64,24) 


1.17 ±0.02 


1.77 ±0.18 


1.76 ±0.19 


Simulation (N= 128,32, 16) 


1.27 ±0.15 


1.72 ±0.10 


1.76 ±0.18 


Charbonneau et al. (2001) 




1.71 ±0.01 


1.74 ±0.06 



Table 3. Theoretically predicted and numerically simulated correlations between the length scale L, time duration T, peak energy 
dissipation rate P, and total energy E for cellular automaton models with Euclidean dimension S - 1 , 2, 3 in the state of self- 
organized criticality. 



Reference 


5 = 1 


5=2 


5=3 




Fractal Dimension: 








Diffusive scaling of time with length: L oc T^''''- = T''^ 








Theory: ySji 


0.50 


0.50 


0.50 


Slopes (N=256,64,24) 




0.68 ± 0.25 


0.49 ± 0.21 


Slopes (N= 128,32, 16) 




0.76 ±0.13 


0.48 ± 0.21 


Regression (N=256,64,24) 


0.53 ±0.15 


0.65 ± 0.02 


0.53 ± 0.02 


Regression (N= 128,32, 16) 


0.65 ±0.13 


0.61 ± 0.01 


0.53 ± 0.02 


Correlation of power with duration: P oc T^'^'' 








Theory: firp 


0.50 


1.00 


1.50 


Slopes (N=256,64,24) 




0.91 ±0.19 


0.79 ± 0.24 


Slopes (N=128,32,16) 




0.99 ±0.12 


0.80 ± 0.22 


Regression (N=256,64,24) 


0.67 ± 0.20 


1.04 ± 0.11 


1.01 ± 0.12 


Regression (N=128,32,16) 


0.73 ±0.14 


1.01 ±0.08 


1.02 ±0.14 


Correlation of energy with duration: E oc T^te 








Theory: firs 


1.50 


1.75 


2.00 


Slopes (N=256,64,24) 




1.61 ±0.08 


1.51 ±0.20 


Slopes (N=128,32,16) 




1.73 ±0.10 


1.48 ±0.19 


Regression (N=256,64,24) 


1.56 ±0.09 


1.80 ± 0.05 


1.75 ± 0.06 


Regression (N=128,32,I6) 


1.65 ±0.07 


1.82 ±0.04 


1.75 ±0.06 



flare durations T, and « 1.56 ± 0.04 for the total energies E 
(Fig. 9, left). In addition we re-analyzed the BATSE/CGRO flare 
data set and obtained similar values of ap a; 1.71 ± 0.06 for the 
peak energy dissipation rate P, aj ~ 2.02 ± 0.28 for the flare du- 
rations T, and oe ~ 1.49 ± 0.07 for the total energies E (Fig. 9, 
right). 



Comparing the observed values with the theoretically pre- 
dicted values of our 3-D fmctal-diffusive SOC model (Section 2 
and Table 1), we find the foUowing ranges regarding the pow- 
erlaw slopes of the occurrence frequency distributions: flare du- 
rations ar = 2.02 - 2.29 (predicted a'^^" = 2.00), peak energy 
dissipate rate orp = 1.71-1.73 (predicted or^''" = 1.67), and total 
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Fig. 9. Occurrence frequency distributions of peak count rates P, durations T, and total counts E in solar flares observed with 
HXRBS/SMM during 1980-1989 (left panels ) and with BATSE/CGRO during 1991-2000 (right panels) at energies > 25 keV (left). 



energies as = 1.49 - 1.56 (predicted a'^'" = 1.50), Thus, the- 
ory and observations agree within a few percents, which is even 
better than the agreement of theory with the cellular automa- 
ton simulations in the 3D case (probably due to the limited grid 
size). The powerlaw index of the peak fluxes was found to vary 
with the solar cycle within a range of arp » 1.6 - 1.9 (Crosby et 
al. 1993; Biesecker 1994; Bai 1993; Aschwanden 2010a), which 
could indicate a variable degree of magnetic complexity that is 
manifested either in a time-dependent fractal dimension Ds , or 
threshold Be of SOC events (Aschwanden 2010a,b). 



Comparing the observed flux time profiles /(/) of the hard 
X-ray flux of large solar flares, they fluctuate erratically in a 
similar way as shown for the largest avalanche of cellular au- 
tomaton simulations (Fig. 3, top left panel), as records with high 
time-resolution and high sensitivity show, e.g., observed with 
BATSE/CGRO (Aschwanden et al. 1998). The erraticaUy fluctu- 
ating hard X-ray time profiles are generally interpreted in terms 
of the chromospheric energy dissipation rate of precipitating 
nonthermal electrons produced in magnetic reconnection pro- 
cesses, similar as we defined the instantaneous energy dissipa- 
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tion rate /(f) (Eq. 8) in our cellular automaton model, while the 
associated soft X-ray time profile represents the thermal emis- 
sion of the heated plasma, which is monotonically increasing 
during the impulsive flare phase because it represents the time 
integral of the heating rate according to the Neupert effect, sim- 
ilar as we defined the total dissipated energy e(f) (Eq. 9) in our 
cellular automaton model. The positive index in the correlation 
of the peak energy dissipation rate P with the time duration T in 
our avalanche model, i.e., P oc T^^^ (Eq. 10), thus predicts that 
the hard X-ray peak flux is statistically the higher the longer a 
solar flare lasts (which is also known as "big-flare syndrome")- 

A specific prediction is that in 3-D Euclidean space, regard- 
less of the fractal dimension in the range of D3 = 1 , 3, the dis- 
tribution of flare energies is restricted to a relatively small range 
of a'^" = 1.40 - 1.67, which is consistent with our observations 
of oe = 1.49 - 1.56. The numerical value of the energy power- 
law slope of aE ~ 1 .5 implies that the total energy of all flares 
is heavily weighted by the largest flares, while nanoflares con- 
tain only an insignificant amount of energy, unless the powerlaw 
slope is steeper than a critical value of =2 (Hudson 1991), 
which is an ongoing argument in the controversy of coronal heat- 
ing by nanoflares. In order to make nanoflares dominant, a pow- 
erlaw slope of qe > 2 is needed, which contradicts our model, 
since the maximum value cannot exceed aE,max = 5/3 = 1.67, 
well below the critical limit of = 2. 

Interestingly, some cellular automaton simulations have 
been performed that actually produced significantly steeper 
powerlaw slopes, say in the order of ~ 3.0, which seems to 
contradict our conclusions. One simulation produced a bro- 
ken powerlaw distribution with a slope of ap - 3.5 for the 
smallest events, interpreted as nanoflare regime, while a flat- 
ter slope of ap - 1.8 was found for the larger flares (Vlahos 
et al. 1995), a difference that is attributed to the anisotropic 
next-neighbor interactions applied therein. Another study 
simulated solar flare events as cascades of reconnecting mag- 
netic loops, with the flnding of a powerlaw slope of = 3.0 
for the released energies (Hughes et al. 2003). Interacting 
loops with a length scale L are bisected during a reconnection 
step into two shorter scales L/2, and the released energies are 
defined in terms of the length scale therein, i.e., E oc L, for 
which our model predicts indeed an occurrence frequency 
distribution of N{E) oc N(L) oc L-^ ^ ^ 3 22) for an 
EucUdean dimension S = 3, so it is fully consistent with our 
model if their energy definition is adopted. Discrepancies of 
powerlaw slopes among different studies can indeed often be 
explained in terms of inconsistent definitions of the energy 
quantity. 

5. CONCLUSIONS 

We developed an analytical theory for the statistical distributions 
and correlations of observable parameters of SOC events, which 
includes the avalanche length scale L, the time duration T, the 
peak P and energy dissipation rate F. The basic assumptions of 
our analytical model, which we call the fractal-diffusive SOC 
model, are the following: 

1. Diffusive expansion of SOC avalanches: The radius r(t) or 
spatial length scale L of an avalanche grows with time like 
the average of a diffusive random walk, which predicts a sta- 
tistical correlation L oc T^l^ between the length scale L and 
time duration T of the avalanche. 

2. Fractal Energy Dissipation Rate: The complexity of ran- 
dom next-neighbor interactions in a critical SOC state can 



be characterized approximately with a fractal geometry. The 
volume (or area) of the instantaneous energy dissipation 
rate is assumed to have a fractal dimension Ds. The pre- 
dicted statistical correlations are: F oc p 

E oc T^+Dsl2_ 

3. Mean Fractal Dimension: The mean fractal dimension Ds 
for different EucUdean space dimensions S = 1,2,3 can 
be estimated from the arithmetic mean of the minimum di- 
mension for a propagating avalanche, Ds,mm ~ 1, and the 
maximum (Euclidean) dimension Ds max - S , which yields 
Ds^(l+S)/2. 

4. Occurrence Frequency Distributions: Equal probabiUty of 
avalanches with size L at various spatial locations in a uni- 
form, slowly-driven SOC system predicts a probability dis- 
tribution of NiL) oc . A direct consequence of this as- 
sumption, together with the other assumptions made above, 
yields a powerlaw function A'(x) oc x'"' for the occurrence 
frequency distributions of all parameters, which is the hall- 
mark of a SOC system. The predicted powerlaw indices are: 
ar = (1+ S)/2, ap = l + (S - l)/Ds, ap = 2- 1/5, and 
oe = I + {S - l)l{Ds + 2). Specifically, for applications to 
3-D phenomena, absolute values are predicted for the power- 
law slopes ai-3 and ap — 2, and ap — \ .(fl, while a range 
of aE = 1 .4, 1 .67 is expected for any fractal dimension in 
the range of 1 < D3 < 3. 



We have validated our theory by a detailed comparison with 
a set of SOC simulations using a specific form of a cellular au- 
tomaton avalanche model (connectivity, stability threshold, re- 
distribution rule, etc), and found a good agreement between the- 
ory and numerical simulations, in the order of ^ 10% for the 
powerlaw slopes (ax, orj, «£,«/>), the power indices of corre- 
lated parameters (fiTL,pTP,PTE), and the fractal dimensions Ds, 
for all three Euclidean space dimensions S - 1, 2, 3. Yet, at its 
most general level our theory is saying that the self-similarity 
of energy release statistics in such models is a direct reflection 
of the fractal nature of avalanches. The SOC flare model re- 
cently proposed by Morales & Charbonneau (2008, 2009) offers 
an interesting test of this conclusion. Their model, defined on 
a set of initially parallel magnetic flux strands, contained in a 
plane and subjected to random sideways deformation, with in- 
stability and readjustment occurring when the crossing angle of 
two flux strands exceeds some threshold angle. This model is 
thus strongly anisotropic, with pseudo-local stability and redis- 
tribution, in the sense that these operators now act on nearest- 
neighbors nodes located along each flux strand, rather than in 
the immediate spatial vicinity of the unstable sites. This is very 
different from the isotropic Lu et al. (1993)-type SOC model 
used here for validation. Yet, the results compiled in Table 2 
of Morales & Charbonneau (2008) for their highest resolution 
simulations reveal that the theoretical occurrence frequency dis- 
tributions obtained herein, do hold within the stated uncertain- 
ties on the power-law fits. Likewise, the power indices of the 
correlated parameters also agree with theory within the inferred 
uncertainties. This provides additional empirical support to our 
conjecture that the fractal-diffusive SOC model does represent a 
robust characterization of avalanche energy release in SOC sys- 
tems in general. However, it should be remembered that the in- 
ferred scahng laws are only valid for a slowly-driven SOC sys- 
tem, while alternative SOC systems with time-variable drivers 
or non- stationary input rates exhibit modified occurrence fre- 
quency and waiting time distributions (Charbonneau et al. 2001; 
Norman et al. 2001), which was also found in solar observa- 
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tions extending over multiple solar cycles (Crosby et al. 1993; 
Biesecker 1994; Bai 1993; Aschwanden 2010a). 

What other predictions can be made from our analytical 
model? For SOC processes in 3-D space, which is probably 
the most common application in the real world, the mean frac- 
tal dimension is predicted to be D3 =s 2.0, which can be tested 
by measurements of fractal dimensions in observations. The 20 
largest solar flares observed with TRACE have been analyzed in 
this respect and an area fractal dimension of D2 = 1.89 ± 0.05 
was found at the flare peaks, which translates into a value of 

- 2.10 + 0.14 if we use an anisotropic flare arcade model 
(Aschwanden and Aschwanden 2008a). The distribution of flare 
energies is predicted to have a powerlaw slope of oe = 1.50, 
which closely matches the observed statistics of solar flare hard 
X-ray emission (as ~ 1.49 - 1.56). Since this value is undis- 
putably below the critical limit a = 2 of the energy integral, the 
total released energy is contained in the largest flares and thus 
rules out any significant nanoflare heating of the solar corona. 
Another prediction, that we did not test here with solar flare data, 
is the diffusive flare size scaling. Straightforward tests could be 
carried out by gathering statistics of the flare size evolution dur- 
ing individual flares (which are predicted to scale as x(t) oc f'^^), 
as well as from the statistics of a large sample of flares, which is 
predicted to show a correlation L oc r'^^. The application of our 
fractal-diffusive SOC model to solar flares imphes that the sub- 
sequent triggering of local magnetic reconnection events during 
a flare occurs as a diffusive random walk. A similar finding of 
diffusive random walk was also found in the turbulent flows of 
magnetic bright points in the lanes between photospheric granu- 
lar convection cells (Lawrence et al. 2001). The spatio-temporal 
scaling of the difliisive random walk predicts also the size, du- 
ration, and energy of the largest flare, which is likely to be con- 

1 /2 

strained by the size Lar « T^^x of the largest active region. 
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